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We report studies of the equilibrium and the dynamics of a general set of lattice models which 
capture the essence of the force-induced or mechanical DNA unzipping transition. Besides yielding 
the whole equilibrium phase diagram in the force vs temperature plane, which reveals the presence 
of an interesting re-entrant unzipping transition for low T, these models enable us to characterize 
the dynamics of the process starting from a non-equilibrium initial condition. The thermal melting 
of the DNA strands displays a model dependent time evolution. On the contrary, our results suggest 
that the dynamical mechanism for the unzipping by force is very robust and the scaling behaviour 
does not depend on the details of the description we adopt. 



The replication of DNA is a correlated process involv- 
ing many proteins and other molecules j^] working at dif- 
ferent points in space and time. An understanding of the 
nature and origin of this correlation is expected to shed 
light on this complex mechanism. It has recently been 
shown (|-|^ that the force induced unzipping of DNA is a 
genuine phase transition different from the thermal melt- 
ing transition of DNA. It was then hypothesized Q that 
the initiation of replication at the origins along the DNA, 
e.g, by dnaA for E.Coli [Qj^ or bythe "origin recognition 
complex" (ORG) in cukaryotcs ||^ is like this unzipping 
near the critical threshold (with dnaA or ORG acting 
as the force-inducing agent) and the resulting correla- 
tion during unzipping leads the co-operativity required 
for replication. 

In contrast to real biological situations, techniques 
like laser tweezers atomic force microscopes (AFM) 
^ |l2| etc have been used to study DNA by pulling at 
one end. This has led to strand separation by force. In 
particular, AFM experiments reported hysteresis in the 
unzipping process, indicating the presence of a first or- 
der transition. These mechanical unzipping experiments 
have opened up new ways of thinking about DNA, just as 
similar stretching experiments of DNA showed the pos- 
sibility of several structures other than the most preva- 
lent B-DNA The activities of polymerases, topoiso- 
merase etc on single stranded DNA have now been ana- 
lyzed in terms of the force they exert or the force applied 
against them [^^-^. What needs to be investigated, to 
mimic the biological situation, is the coupling between 
the opening of the strands and the subsequent events 
during replication. Such a study involves the dynamics 
of the unzipping process |^ . 

The purpose of this paper is to define a set of sim- 
pler models, in the spirit of Poland and Sheraga jl7|, for 
which the unzipping transition can be studied exactly. 
On the basis of this, the dynamics can be understood. 
The proposed lattice models (bubble models: b-models) 



incorporate the mutual-avoidance (hard-core repulsion) 
of the strands (and also self-avoidance). A further simpli- 
fication is obtained by suppressing the bubbles along the 
chains, thereby defining a "fork model" or "Y-model". 
The phase diagram of the equilibrium system displays 
in both cases a re-entrant region at low T: for a finite 
range of forces the molecule gets unzipped by decreas- 
ing the temperature. The dynamics of both the b- and 
the Y-models in the various phases and on the phase 
boundary are then studied, by starting from a zig-zag 
non-equilibrium bound state as the initial condition. We 
find that in all regimes above or below the critical line 
in the phase diagram, the time evolutions of the order 
parameters follow dynamical scaling laws. 



The models and their equilibrium phase diagram: 

We model the two strands of DNA by two directed self 
and mutually avoiding walks. In two dimensions, on the 
square lattice (see Fig.l), this means that the two walks 
are forced to follow the positive direction of the diagonal 
axis (1,1) (i.e. the coordinate along the direction (1,1), 
to be called the z direction, always increases). The force 
is then acting along the direction perpendicular to direc- 
tion z, i.e., (—1, 1): this transverse direction is called the 
x-direction. By measuring the x separation in unit of the 
elementary square diagonal, we say that two complemen- 
tary monomers are in contact when x = 1: a (positive) 
binding energy e is gained for each contact. Due to the 
geometrical properties of the lattice, all these contacts 
contributing to the energy involve monomers labelled by 
the same z-coordinate, as one would require for base pair- 
ing in DNA (see Fig. 1). 
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and it has been plotted in Fig. 2. 




FIG. 1. A typical configuration of tlie two DNA-strands as 
modelled on a square lattice. Dashed lines indicate monomers 
which are in contact. The quantities m and x are graphically 
represented. 

By fixing e = +1 throughout the calculations, all the 
thermodynamic properties of the system depend on the 
temperature T — (3~^ and on the force /. We have set 
fcs = 1. The partition function of the two n-step chains, 
one starting at (0, 1) and the other at (1, 0) is 



exp {(3 fx) 
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where dn (x) represents the fixed distance partition func- 
tion, i.e., the sum over all interacting pairs of directed 
chains whose last monomers are at distance x fl^ and 
T is the free energy density. The quantity dn{x) whose 
/3 dependence has been omitted, obeys simple recursion 
relations. In d = 2, on the square lattice, the recursion 
relation is 

dn+l{x) = [2dn{x) + dn{x + 1) + dn{x - 1)] 

X [1 + (e^ - 1) , (2) 

valid for x > 1 with the conditions d^^x) = S^^i and 
dn{0) — for all n. Note that this kind of equation 
appears also in other models of DNA [|l7|jl9[ and in stud- 
ies of random walks adsorption The model can be 
asymptotically solved by locating the singularity closest 
to the origin of its related generating function 



(3) 
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A phase transition is indicated by any singular change in 
the location of the closest singularity. 

With some simple calculations a critical line, fc{T), 
separating the zipped from the unzipped phase can be 
obtained 



/,(T) = Tcosh- 
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FIG. 2. The plot of the force vs. temperature phase dia- 
gram for the model on the directed lattice. 



At / = the critical melting temperature is Tm = 
log 4/3 ' critical force at T = is /o = 1 whereas 
for / > /max = 1.358806... the system is always in the 
unzipped state. Similar results can be obtained for the 
physically relevant three dimensional case and for higher 
dimensions |^ . The average fraction of contacts Q is zero 
(not zero), when / > fc{T) (/ < fc{T)) while the average 
end-to-end distance {x) = —N^, where N is the num- 
ber of base pairs in each chain, and the average number of 
"liberated" monomers {i.e. from the last contact to the 
end), m (See Fig. P have the opposite behavior. The 
transition line is first order at any point with fc(T) ^ 
in d = 2, 3, 4, 5 and everywhere for d > 5. 
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We notice that there is a re-entrance in the / — T phase 
diagram. For /max > / > /o the usual denaturation tran- 
sition is present but if the temperature is further lowered 
the two strands separate again through a "cold unzip- 
ping" 11,11 . 

A further simplification to this model can be obtained 
by suppressing the bubbles along the chain, i.e. by 
considering only conformations having the first N — m 
monomers bounded, whereas the remaining m are sepa- 
rated in a Y-like conformation. 

This Y-model will be extremely useful to study the un- 
zipping dynamics and it presents a phase diagram simi- 
lar to the one previously obtained but, for example, with 
T,n = Ts|2' /o = 1 and /max = 1.282143... [|4) in d=2. 

Dynamics. We now consider the dynamics of the mod- 
els previously introduced, the b-model of directed walks 
and the simplified Y-model. In both cases, we start from 
a non-equilibrium "zero-temperature" initial condition 
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with the two chains zipped in a zig-zag configuration, 
and let the system evolve at a temperature T and under 
a force /, with T and / chosen so that the equilibrium 
state is either on or above the critical line. The simpler 
2-dimensional case will be presented first and the gener- 
alization to higher dimensions will be discussed later. 

The five regimes considered are marked A-E in Fig. 
Ij. Numerically, a Monte-Carlo dynamics is used to mon- 
itor the time evolution of the order parameters m and 
X previously introduced (see also Fig. 1). For the dy- 
namics, we selected only local physical moves, so that 
the model should be the lattice counterpart of the Rouse 
model in the continuum. In all cases we find that, far 
from saturation, the order parameters evolve as power 
laws as functions of time t: in particular we can define 
two exponents 6i and 62 as 



m(f) 
JV 



for various values of A^: the col- 



i{t) - t^^ , and x{t) - .t(0) ^ 
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Notice that m and x in Eq. ^ refer to the ensemble 
average values, but average signs have been omitted for 
simplicity of notation. We in fact find the following dy- 
namical scaling laws to hold: 



m(t) 



t \ x{t) - a;(0) 
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where N is the length of the chains and Gm.x are two scal- 
ing functions. Eq. ^ also defines the exponents di_2 and 

zi 2 for the two variables. Note that di ~ —, 62 = — and 

. . . ^1. -^2. 

that di 2 can be obtained through equilibrium consider- 
ations as one requires m,{t) iV^^ and x{t) N'^^. 
The crossover to the equilibrium behavior is described 
by the "dynamic" exponents zi.2. The possibility of the 
two quantities having different relaxations is kept open 
by two different exponents. 

The values of the exponents in Eqs. (|^) and ^ ob- 
tained from simulations, for both the Y-model and the 
b-model, are shown in Table I. They were obtained by 
collapsing Monte-Carlo data according to e q. (m ) and by 
using a recently proposed search algorithm EH] (see Fig. 
3) 
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lapse of all curves indicates that d\ = 1, z\ 



§ and 61 = |. 



It is possible to explain the exponents for the Y-model 
found numerically. There are two relevant mechanisms 
which can drive the separation of the strands: the un- 
binding of bases at the bifurcation point and the stretch- 
ing of the strands at the extremes. The combination of 
these two processes controls the overall behavior in dif- 
ferent regimes. 

Regime A: f = 0, T > Tm- Above the critical tem- 
perature the dominance of the entropy implies that at 
every time step one base pair breaks, yielding a linear 
behavior with 9i = 1 and di = 1. Also x{t) tends to in- 
crease, up to its equilibrium value A'' 2 : this is reflected in 
the equilibrium probability distribution P^q (x) which dis- 
plays an upward derivative aX x = 1. This suggests (see 
p5[ ) that the dynamics of this quantity should be in the 
same universality class of the one-dimensional Kardar- 
Parisi-Zhang equation p5| , and so ^2 = 3, 0^2 = 5, and 

Z2 = |. 

Regime B: f = 0, T = T„j: In this regime ^1 = 5 = 
z^^ because at criticality the probabilities to increase 
and to decrease m are expected to be equal, so that m 
performs, roughly speaking, a random walk in time with 
reflecting boundaries at m = 0, N. Also for the end-to- 
end distance, steps toward larger or smaller values of x 
are equally probable, and this means that the equation 
obeyed by x should stay in the same universality class 
as the one-dimensional Edward- Wilkinson equation, and 
therefore O2 = j and Z2 = 2. 
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Regime C: f > 0,T > Tm- The strands microscopi- 
cally tend to stretch along the direction of the pulling 
force. However, once we have pulled the two chains up 
to an end-to-end distance x, to increase x further by one 
unit we flrst need to move all the stretched part, which 
would take a time typically of order x. In other words, 
one has x{t + tg) ^ x{t) + ^ and we expect the dynam- 
ical exponent 62 to be ^ (of course 6'i = 1 as before and 
di=d2 = 1). 

Regime D: T < Tm,f > fc{T): Here the only micro- 
scopic mechanism for opening the fork is through the 
applied force: the strands must stretch completely in the 
vicinity of the bifurcation point and only at this point 
will the fork liberate one more monomer, because oth- 
erwise the opening is energetically very unlikely. Thus 
we expect that x m and, using arguments as done for 
regime C, 0i = 6*2 = ^, di = ^2 = 1- In ref. |^ it was 
found, in the mean fleld approach and in a model resem- 
bling our Y-model, that the time necessary to unzip the 
two strands completely is of the order of A^. This is 
in agreement with our analysis, but works only in this 
regime. 

Regime E: T < T^J — fc{T): At criticality, one 
expects that the cost for unzipping and zipping is the 
same (the equilibrium probability distribution of having 
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m monomers unzipped or an end-to-end distance equal 
to X is flat), therefore x{t + to) = x{t) ± ^ with equal 
probability. Therefore, the end-to-end distance makes a 
random walk in the rescaled time j so that x ~ (j)^ 
and 02 (= 01 ) = 5 and di = (i2 = 1 since at coexistence 
there is a finite fraction of liberated monomers. 

Another way of obtaining 0i = i is to demand that a 
kink liberated at the fork needs to diffuse out of the end 
before the next one is released. In other words, the rate 
of change of m is determined by the diffusion of a kink 
over a distance m. The latter time-scale being of order 
TO~^, we expect dra/dt ^ which gives m ^ . 

This simple model deserves some observations. Re- 
markably, it displays two different time scales in its dy- 
namics, mirrored in the difference of the exponents zi^2 
in regimes A and C. The first time scale {t ^ N^'^) quan- 
tifies the time necessary for the unbinding (or unzipping) 
of the bases while the second {t ^ N^^) gives the time 
needed to open (and to stretch whenever f 0) the 
two chains up to their equilibrium end-to-end distance. 
At T < Tjn the two processes are virtually the same, 
because the unbinding (or unzipping) is dragged by the 
stretching. However, above T^, the processes decouple 
and the unbinding gets faster, being controlled by the 
temperature, and therefore zi < Z2. 

Furthermore, in the numerical calculations we found 
large sample-to-sample fluctuations, thereby requiring a 
huge number of runs to get good averages of m{t) and 
x(t). This is due to the long time correlation that exists 
in the system, which keeps samples with different initial 
histories far apart for any t. 

Turning to the b- model with bubbles, amazingly the 
dynamical exponents in regimes B, E and D with T < T„i 
are the same as found in the Y-model. This establishes 
that at T < Tm not only for statics, as we saw in the pre- 
vious section, but also for the dynamics, bubbles are not 
relevant in the scaling properties (for T = Tm the equal- 
ity is due to another reason ||2^). At T > T„i, on the 
other hand, the opening of bubbles heavily affects the 
base unpairing process, unlike the Y-model case where 
bubbles are forbidden. The length of the unzipped part 
in the present case now can change by ±Z(i), where l(t) 
is the typical length of bubbles, and the motion of the 
fork point can by no means be approximated by a simple 
random walk ( and so 9i changes as shown in Table I). 
The quantity x{t) instead has a dynamics in the b- model 
similar to the fork case and indeed 02 is the same for both 
the models in all regimes. We show in Fig. 3 the collapse 
leading to 9i — j in regime A. 

An important question is that of the dependence of 
these results on dimensionality. As for the Y-model, 
the arguments we gave suggest that there be no d- 
dependence. For the b-model, instead, this should be 
true only in regimes D and E for T < Tm, where the 
Y-model gives the exact result; at T > Tm, on the other 
hand, bubbles play a dominant role and so we expect a 
dependence on d. We confirmed this picture with some 



calculations on a simpler model which should be in the 
same universality class of the one under study: that of a 
single random walk, pinned at the origin by an attractive 
interaction and subject to a stretching external force. In 
this system m is defined as the number of monomers from 
the last visit to the origin to the end of the walk. For 
T > Tm, in regimes A and C, our calculations show that 
the exponent 6i increases as dimension increases, appar- 
ently with no upper critical dimension. Just at criticality 
at zero force, instead, we find that the exponent 6i is very 
close to ^ in any d. The emerging picture of robust results 
for T < Tm and model-dependent dynamics for T > Tm 
would be preserved even if, in the original models, the 
directedness constraint is relaxed (as is the case for the 
statics p3|). 

The arguments presented so far could be generalized 
to add other ingredients of dynamics also. An example is 
the nonlocal effects in dynamics, as is the case e.g. in the 
Zimm model. In ref. ^J, the author suggests that, in the 
regime we call D, in the mean field approximation, non- 
locality could speed up the opening so that 0i^2 would 
be instead of ^ (m is the exponent characterizing the 
length dependence of the mobility as defined, e.g., for the 
Zimm model, see js) for the notation and [g9| for a re- 
view). In regime E, we can combine our arguments with 
the same reasoning to get 9i ~ 62 = As we expect, 
the Rouse model results are obtained with u = 1. 

Lastly, the above analysis can also be extended to 
binary disorder in the contact potential, which is of 
course a realistic feature to be included in the model, 
i.e. the energy of the contact of the z-th base may 
be = e + Ai, being a random variable with bi- 
nary distribution and zero mean. The Y-model offers 
a good starting point for the study of the effects of 
heterogeneity, because the critical line of the quenched 
model can be proved to be the same as that of the pure 
model with energy e. If we call !F(m) the free energy 
density of a configuration with m (out of N) unzipped 
monomers, we expect on general grounds that P{m,t), 
the probability of having m monomers unzipped at time 
t, will obey a master equation with transition rates W±^m 
depending on the realization of disorder, i.e. we ex- 
pect that W±^m oc minje^F'^^^^^,!!, with NP^^^ 

= + log(l -I- cosh(|)) II, which can be 

seen to consist of a zero- mean random "noise" ( ) 
plus a "bias" . We thus believe that the dynamics of the 
Y-model can be mapped onto Temkin's model of a ran- 
dom walker in a random environment pO| (provided that 
at r < Tm we rescale time as described above in regimes 
D and E ), so that, following Ref. |^^, there will be 
a region around the critical line, at T < Tm, in which 
m,x t2 (and just on the critical line at / 7^ one 

is tempted to expect m,x ^ (log (j)) ). The quantities 
m, X are now quenched averages over realizations of dis- 
order. The scenario we propose is sketched in Fig. 4. The 
curve bounding the region where, at T < Tm, disorder 
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should be important for the dynamics has been found 
by applying the criterion of Ref. |^ to the system with 
the above transition rates. This curve coincides with the 
critical line of the static model with annealed disordered. 
Preliminary runs at very low T suggest that as T ^ 0, 
for a given realization of disorder and with any force be- 
tween ei and £2 (the two binding energies), the unzipping 
will take place as far as the first more attractive base 
pair is found. Besides, the probability that the monomer 
of index toq will ever be unzipped decays exponentially: 
P(mo) c!c 2"™°. This is consistent with the mapping sug- 
gested above. A detailed study of the effect of disorder 
will be done elsewhere. 

Conclusion. We obtained the phase diagram of a lattice 
model for the unzipping transition of a double stranded 
DNA by a force. This model incorporates both the mu- 
tual and self-avoiding nature of the two strands and the 
equilibrium problem can be solved exactly. A still simpler 
model has been defined by suppressing the bubbles that 
are important for thermal denaturation. This simpler 
version is also exactly solvable and is shown to retain the 
basic features of the b-model. Thanks to the exact knowl- 
edge of the phase boundary, we have been able to inves- 
tigate the dynamics of unzipping from a non-equilibrium 
bound state both on and away from the phase bound- 
ary. The dynamics shows scaling behaviors in different 
regimes of the phase diagram. These scalings in most 
cases could be understood from the plausible mechanisms 
of unzipping and denaturation as discussed in the previ- 
ous sections. We end with two notable features. In the 
case of the b-model with bubbles, the unzipping at high 
temperatures remains a puzzling issue, especially the dy- 
namic exponent zi = | with or without force. We believe 
that this has to do with the statistics of bubbles, though 
no satisfactory answer could be found. This remains an 
open problem. Lastly, on a positive note, the unzipping 
dynamics on the phase boundary in the presence of a 
force is found to be distinctly different from the thermal 
denaturation at zero force. Whether a real biological sys- 
tem takes advantage of these differences to distinguish 
the unzipped region of DNA from a fluctuation-induced 
bubble formation remains to be probed. 



FIG. 4. Plots of the phase diagram for homo-DNA (top) 
made up of the more attractive base (e2 > ei), and for the 
annealed (middle) and quenched (bottom) disordered models. 
The numbers refer to the values of the exponent 6\ in the var- 
ious regions of the phase diagram, which we believe to hold 
on the basis of the mapping onto Temkin's model. Note that 
above the annealed critical line at T < and everywhere at 
T > Ta the pure system results are recovered. Interestingly, 
the annealed line does not show the re-entrant behavior found 
in Section II and in the quenched system. "?" indicates the 
regime where we do not have numerical evidence in support of 
our results. We took t2 = 1.5 and ei = 0.5 in the calculations. 
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TABLE I. "Dynamic" and equilibrium exponents for the 
Y-model (Y) and the b-model with bubbles (b) as defined in 
eq.(0). The regimes A,B,C,D and E are those shown in Fig.l. 
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